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Abstract: No attempt has been made to date to model growth in girth of 
rubber tree (Hevea brasiliansis). We evaluated the few widely used 
growth functions to identify the most parsimonious and biologically 
reasonable model for describing the girth growth of young rubber trees 
based on an incomplete set of young age measurements. Monthly data for 
girth of immature trees (age 2 to 12 years) from two locations were sub¬ 
jected to modelling. Re-parameterized, unconstrained and constrained 
growth functions of Richards (RM), Gompertz (GM) and the monomo- 
lecular model (MM) were fitted to data. Duration of growth was the 
constraint introduced. In the first stage, we attempted a population aver¬ 
age (PA) model to capture the trend in growth. The best PA model was 
fitted as a subject specific (SS) model. We used appropriate error vari¬ 
ance-covariance structure to account for correlation due to repeated 
measurements over time. Unconstrained functions underestimated the 
asymptotic maximum that did not reflect the carrying capacity of the 
locations. Underestimations were attributed to the partial set of meas¬ 
urements made during the early growth phase of the trees. MM proved 
superior to RM and GM. In the random coefficient models, both Gf and 
Go appeared to be influenced by tree level effects. Inclusion of diagonal 
definite positive matrix removed the correlation between random effects. 
The results were similar at both locations. In the overall assessment MM 
appeared as the candidate model for studying the girth-age relationships 
in Hevea trees. Based on the fitted model we conclude that, in Hevea 
trees, growth rate is maintained at maximum value at to, then decreases 
until the final state at dG/dt > 0, resulting in yield curve with no period of 
accelerating growth. One physiological explanation is that photo synthetic 
activity in Hevea trees decreases as girth increases and constructive 
metabolism is larger than destructive metabolism. 
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Introduction 

Hevea brasiliensis (Willd. ex. Adr. de. Juss.) Muell. Arg., is an 
important industrial tree crop grown mainly in tropical climates 
up to 12° latitude on both sides of the equator. In India, the crop 
is largely cultivated between 8° 15' N and 12° 52' N latitudes. 
Hevea rubber is the major source of ‘natural rubber’, a product 
used in the manufacture of thousands of products of which the 
pneumatic tyre is the most common. Rubber trees belong to the 
genus Hevea of the family Euphorbiaceae. Though ten species 
have been recognized in the genus, H. brasiliensis is the only 
species cultivated for obtaining rubber, a product of commercial 
importance (Webster and Paardekooper 1989). The trees are 
erect with straight trunks, growing in the wild to over 40 m with 
a life span of more than 100 years. In managed plantations they 
rarely exceed 25 to 30 m because tapping reduces growth and 
trees are usually replanted after 25^0 years when yields fall to 
an uneconomic level. Most of the modern plantations are raised 
by clonal buddings on seedling rootstocks. Normally the trees are 
tapped for latex after 5-8 years of growth and tapping continues 
for about 204>5 years after which the trees are cut and new plan¬ 
tations are raised. Rubber trees are immature until 6-8 years of 
age. During this stage, trunk girth measurements and the calcu¬ 
lated girth increments are widely used as parameters of growth. 
The trunk girth and rate of increase in girth are used in experi¬ 
mental work to assess the growth performance of new planting 
materials and the effects of cultural treatments on growth (Shor- 
rocks et al. 1965; Chandrashekar et al. 1998; Chandrasekhar et al. 
2005; Chandrasekhar 2007). Similar to other species, increased 
girth in Hevea can be related directly to the accumulation of dry 
matter, which reflects the fraction of growth remaining after 
respiration. Even though domestication of Hevea took place 
more than 150 years ago and Hevea is now grown in diverse 
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environments (Baulkwill 1989), no attempt has been made to fit 
different growth functions and recommend a suitable model ap¬ 
propriate for the crop. 

Biological growth, which is regulated by genetics and envi¬ 
ronment, changes continuously with age, leading to an increase 
in volume, size or shape of organism. Several types of age- 
growth models are widely used to simulate the development of 
individual trees and to model the dynamics of growth (Rat- 
kowsky 1988; Philip 1994; Vanclay 1994; Seber and Wild 2003). 
The standard way to study these changes is to fit mathematical 
functions to growth measurements collected at regular intervals. 
These functions are explicit expressions of both theory and em¬ 
pirical knowledge and are useful for growth and yield prediction, 
computation of sustainable yield, and evaluation of alternative 
management strategies among other uses (El-Shaarawi and Pie- 
gorsch 2002). When growth is described by an appropriate func¬ 
tion, this can help condense information contained in long data 
series into a few parameters from which biologically meaningful 
constants can be derived. From the fitted models, parameters of 
growth can be derived, including maximum growth potential, 
maximum rate of change, and site index of clones (a parameter to 
indicate performance of a clone on a site for a chosen reference 
age),. These parameters are of importance for assessment of the 
performance of clones across environments. It is most desirable 
if these parameters are derived from a set of measurements made 
during the early growth phase without sacrificing accuracy rather 
than from a set of size-age measurements from many trees at all 
ages until the end of the physiologically active life span of the 
species. The objective of the present study was to evaluate a few 
widely used growth functions and to identify the most parsimo¬ 
nious and biologically reasonable model that best describes the 
relationship between girth and age in rubber trees given an in¬ 
complete set of young age measurements. 

Materials and methods 

Study locations, plant material and experimental details 

Data for this study were collected from two contrasting locations 
(LI and L2). LI was the headquarters of the Rubber Research 
Institute of India at Kottayam in Kerala State. L2 was the Rubber 
Research Institute of India’s Regional Research Station at Dap- 


chari in the Thane District of Maharashtra State. The experimen¬ 
tal site at LI was on a hill slope, with Ustic Kanhaplohumults 
soil of bulk density 1.25 g-cm" 3 , permanent wilting point at 
19.7% and field capacity at 27.8%. At L2, data were collected 
from a clone evaluation trial whose site was slightly undulating. 
Soil at this site was a clay loam type with bulk density of 1.30 
g-cm' 3 , permanent wilting point at 17.5% and field capacity at 
30.5%. 

The clone used in the study was RRII 105, the clone most 
widely planted in India. It is a secondary clone produced in India 
by crossing Tjir 1 and G1 1 that are primary clones from Indone¬ 
sia and Malaysia, respectively. Data from LI were from a trial 
where 13 modem clones were under evaluation while at L2 15 
modern clones were under evaluation. At both locations the trials 
were laid out using randomized complete block design. The trial 
at LI had seven replications with experimental plots of seven 
plants in a row in contour planting at a spacing of 3.4 m x 6.7 m 
equivalent to a plant density of 445 plants-ha" 1 . The trial at L2 
had three replications with experimental plots of 36 plants in 
square planting at a spacing of 4 m x 4 m giving a plant density 
of 638 plants-ha' 1 . Geographical, weather, and soil characteristics 
of the study locations are given in Table 1. At LI, field planting 
of the trials was carried out during July/August 1989 while at L2 
it was done during 1985. Three whorled plants raised in polyeth¬ 
ylene bags (55 cm x 25 cm, lay flat dimension) for five to six 
months were used for planting. Standard cultural practices, in¬ 
cluding manuring, weeding, mulching, and sun scorch prevention, 
were all part of agromanagement (Potty et al. 1980; Pushpadas 
and Ahmmed 1980; Punnose et al. 2000). 

Growth data 

Girth measurements made at a fixed height of 1.5 m (the stan¬ 
dard in Hevea research) from the bud union (collar region) were 
subjected to modelling. In each month, the measurements were 
taken during the first three days and were considered as repre¬ 
senting the girth of the trees at the end of the previous month. 
Complete details of the collected growth data are provided in 
Table 2. At LI data were collected from February 1989 through 
January 1997 while at L2 data were collected from February 
1991 through January 1997. A total of six years of data from LI 
and eight years of data from L2 were used in this study. 


Table 1. Geographical, weather and soil characteristics of the study locations 


Location 

Geographical 

Longitude 

Altitude 

Cultivation zone 

Mean annual 

Wet period 

Dry period 

Annual range of 

Latitude (°N) 

(°E) 

(m, asl) 


rainfall (mm) 

(months) 

(months) 

sunshine duration (h) 

Location 1 

9.32 

76.36 

73 

Traditional 


April to Nov. 

Dec. to March 

1.7 to 10.2 

Location 2 

20.04 

72.04 

48 

Non-traditional 


June to Sep. 

October to May 

3.9 to 10.4 


Mean annual 

Mean annual 

Mean annual 


Permanent 

Field capacity 

Bulk density 


Location 

evaporation 

minimum tempera- 

maximum tern- 

Soil Type 

wilting point 

(%) 

(g-cm' 3 ) 



(mm/day) 

ture (°C) 

perature (°C) 


(%) 




Location 1 

3.7 

23 

31.6 

Laterite 

19.7 

27.8 

1.25 


Location 2 

4.5 

20.4 

32.7 

Clay loam 

17.5 

30.5 

1.3 
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Table 2. Details of data for analysis from the two locations 



Location 1 

Location 2 



Number of 






years for 



Number of years 



which data 



for which data is 

Year 

Remarks Age is available 

Remarks 

Age 

available 




Year of 



1985 



planting 

1 


1986 




2 


1987 




3 


1988 




4 


1989 

Year of 


Start of data 




planting 

1 

collection 

5 

1 

1990 


2 

— 

6 

2 

1991 

Start of 






data 






collection 

3 1 

— 

7 

3 

1992 

— 

4 2 

— 

8 

4 

1993 

— 

5 3 

— 

9 

5 

1994 

— 

6 4 

— 

10 

6 

1995 

— 

7 5 

— 

11 

7 

1996 

— 

8 6 

— 

12 

8 

1997 

•January 

was the end point 

•January was 

the end point of girth 


of girth recording 

recording 




• Tapping started from March 

•Tapping started from October 


Growth functions used 

Many mathematical models are used to describe the growth of 
plants and animals over time. The models are based on general 
deterministic differential equation of the form dy/dt=f(y,t). 
Growth functions used in this study are the Richards, Gompertz 
and monomolecular. These functions belong to the exponential 
decline family of growth models that seek biological interpreta¬ 
tion. Basic constructs of the functions and their mathematical 
properties are given in Table 3. In all three models G is the ob¬ 
served girth at age t (expressed in this study in years), G f is the 


asymptotic limit when t approaches infinity and it indicates the 
maximum possible value of the dependent variable determined 
by the productive capacity of the site, ‘m’ in the Richards model 
is the dimensionless shape parameter that permits a variable 
point of inflection, “k” is a function of the ratio of maximum 
growth rate to mature size, normally referred to as maturing rate. 
It serves both as a measure of growth rate and of rate of change 
in growth rate. Parameter “b” has no specific biological meaning 
and it positions the curve in relation to the time axis (origin of 
the growth curve). In this study, models were reparameterized in 
terms of G 0 = G(t 0 ) for some t 0 were used (Seber and Wild 2003; 
Thornley and France 2007). The models were: 

Richards (RM 0 ) as: 

i 

G = (Gf 1_m) - (G|. 1_m) - Go” m) )e“ k(t ~ to) ) (1_m) 


Gompertz (GM y ) as: G = G f e 


In 



e - k(t - t o ) 


Monomolecular (MMu) as: G = (G f — (G f — G 0 )e to ^) 


These reparameterized models were preferred over basic 
forms because specific meanings could be assigned to each of the 
parameters. The Gompertz equation assumes that the substrate is 
non-limiting and quantity of growth machinery is proportional to 


G with inflection point fixed at 



. Monomolecular model is the 


least complex of the equations that describe the progress of 
growth in a simple, irreversible first order reaction. The model 
possesses only an asymptote and intersects the time axis but it 
has no inflection point. This curve rises rapidly during the initial 
period and the instantaneous growth rate decreases monotoni- 
cally to approach asymptotically some final value. 


Table 3. Basic forms of the functions used in the study and their mathematical properties 


Property 

Equation 

No. of constants 
Growth rate (1 st Derivative) 

Asymptotes - Lower 
Upper 

Curve shape symmetry 
Inflection point 
Maximum growth rate 


Richards (RM) 1 
G = G f (l - be ~ b )i^ 


Functions 


Gompertz (GM)‘ 


G 


G f e 


- be 


- kt 


dG 

dt 


rjG m 


4 


kG 


For m ^ 0 both are asymptotes, but for 

m < 0 there is no lower asymptote and 

no point of inflection 

Asymmetrical 

l 

G f m l ~ m 

m 

G f km l ~ m 


dG 

dt 


3 

G f 

kG log —- 
G 

0 



Asymmetrical 



e 


^Richards, 1959); 2 (Winsor, 1932; Richards, 1959) 


Monomolecular (MM) 1 
G = G , (l - be ~ h ) 


dG_ 

dt 


3 

k(G f - G) 


nil (crosses t-axis) 



Concave 


Absent 


Occurs in the beginning 
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Constraining of functions 


When the models were fitted using unconstrained functions, the 
models underestimated the asymptotic maximum that did not 
reflect the carrying capacity of the locations. Published research 
suggests constrained modelling to circumvent the problem of 
underestimation (Shifley and Brand 1984; Salas and Garcia 
2006). In this study, constraint was put on growth duration (T). 
From previous experience in typical commercial situations, un¬ 
tapped bud-grafted trees take a growth duration (T) of about 25 
to 30 years to attain asymptotic maximum growth (Joseph G. 
Marattukalam, Retired Scientist, Marattukalam House, 225 Ara- 
manapady, Chenganachery-686101, Kerala, India, personal 
communication). A higher range of T= 30 years was considered 
in this study to allow trees to reach their potential size deter¬ 
mined by genetic and environmental factors. The constraint was 
introduced into the growth functions by replacing the rate con¬ 
stant ‘ k ’ in the equations by a rearranged formula that is used 
for working out duration of growth (T) which is a function of 
‘ k ’. This was a simple extension of the ideas that have been 
used to solve for various growth quantities (Richards 1959). The 
formulas were: k=(2m+2)/T, for the Richards model, k=4/T, for 
the Gompertz model and k=2/T, for the Monomolecular model. 
Constructs of the equations used for the re-parameterized con¬ 
strained models were: 

Richards (RM C ) as: 

( 2m+2 Y 1 

G = (GY m) -(GY m) -G[, 1 - m) )Y'Y r / ^ 


Gompertz (GM C ) as: 



G f e 


- In 


f G 

G 

V 


/ 

0 



J J 


Monomolecular (MM C ) as: 

G = (G f -(G f -G 0 )e"T t " t ° ) 


Statistical analyses 

All modelling work was undertaken at the location level using 
the individual tree measurements from all the replicate plots. 
Only tree level structure was considered and the plot level struc¬ 
ture was ignored. Data analyses were carried out in the open 
source software R, a free software environment for statistical 
computing and graphics (R-development core team 2010). A 
two-stage approach was adopted. In the first stage, growth was 
modelled to capture the time trend around a population average 
(PA) model wherein the parameters of the function are assumed 
to be invariant across the trees (i.e., fixed). Fittings of functions 

were carried out using the nlstools package (Delignette- 
Muller and Baty 2010). Starting values for all the parameters 
were specified following Fekedulegn et al. (1999), Ogle (2010). 

Value of t Q was fixed at 7, as this is the standard age for decid¬ 
ing the tappability of Hevea plantations (G 0 at t 0 =7) for latex 
4^ Springer 


harvesting. This also helped in comparison of predicted girth at 
that age in different models. Starting values were first evaluated 
in the function preview of the nlstools package with repeated 
modifications to reach a good approximation of estimates which 
were then used for fitting the model in the nlstools package 
(Delignette-Muller and Baty 2010). Statistical significance of the 
parameters of the nonlinear models was determined by evaluat¬ 
ing the 95% asymptotic confidence intervals (ACI) of the esti¬ 
mated parameters. The null hypothesis (H 0 ) that the parameter in 
question is not significantly different from zero was rejected 
when the 95% ACI of the parameter did not include zero. 

General goodness of fit of the fitted models were evaluated us- 

ing adjusted proportion of variation ( R ac y ), residual sum of 

squares (RSS), Akike information criterion (AIC), Bayesian 
information criterion (BIC) and Model Selection Criterion (MSC, 
Scientist software, Micromath, Saint Louis, Missouri, USA). 
Using AIC and BIC values, the more appropriate model was 
judged as the one with smaller values. The decisions rule with 
MSC is that the model with larger MSC is the more appropriate 
one. Validity of the assumptions of the error model was checked 
by plotting non-transformed residuals against fitted values. Boot¬ 
strap (with 2000 iterations) sampling was also performed to vali¬ 
date the model parameter estimates following Delignette-Muller 
and Baty (2010). 

In the second stage mixed model analysis was attempted to fit 
a subject specific (SS) model (random coefficient regression 
model) allowing for individual differences G f and G 0 . This analy¬ 
sis was restricted only to the PA model selected in the first stage. 
Usually forest growth and yield data exhibit both heteoscedastic- 
ity and autocorrelation (Gregoire et al. 1987). As the data used in 
this study are repeated measurements over time, it is expected 
that measurements on the same tree made near in time would be 
more highly associated than measurements made further apart in 
time. To account for this variation, different error variance- 
covariance structures were tried but only the relevant and helpful 
ones are shown and discussed here as the other methods did not 
yield useful results that warranted presentation. These analyses 
were carried out following Gregoire et al. (1987), Vonesh and 
Chinchilli (1997) and Pinheiro and Bates (2000). 


Results 

Unconstrained/constrained models 

Fitting results of the unconstrained functions (Table 4) revealed 
that parameter estimates for all the functions are significant at a 
= 0.05. Model fittings of LI yielded the lowest estimate for G f in 
GM U? the highest estimate in MM y and an intermediate in RMy. 
The range of 95% ACIs was similar for RMy and MMy but in 
GMu it was small. Estimates of G 0 , ASEs, and 95% ACIs were 
similar for all models. The estimate of m in RMy was a positive 

value with a very high correlation with G f . RMSE, R a( y and 

MSCs were almost similar. L2 yielded the highest value in RM U? 
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lowest in GMy and was intermediate in MMy. Asymptotic error 
and range of 95% ACI was highest in RMy while they were 
similar in the other two models. Estimates of G 0 were similar in 
all models. At L2, m in RMy was negative with a very high cor¬ 
relation with G f . Other fit statistics were almost similar except 
for GMy where AIC and BIC values were higher. The trajectory 
of the fitted curve appeared more realistic in the scatter plots of 
girth of trees overlaid with fitted curves (Fig. 1). Residual plots 
(Fig. 2) did not reveal any non-randomness in errors as they were 
more or less uniformly distributed on either side of mean zero. 

100 
80 
60 
40 
20 
0 

100 
80 

I 60 

§ 40 

20 
0 

100 
80 
60 
40 
20 
0 

Fig. 1 Scatter plot of girth of trees overlaid with function plots of 
unconstrained RM (A, B), GM (C, D) and MM (E, F) models fitted to 
data from two locations LI (A, C, E) and L2 (B, D, F) 

In the results from constrained models (Table 5), the maxi¬ 
mum estimate was obtained in GM C , the minimum in MM C and 
an intermediate in RM C for LI. ASE and range of 95% ACI was 
lowest for MM C while in the other two they were higher. Esti¬ 
mates, ASE’s, and 95% ACI’s for G 0 were almost similar in all 
the models. Estimate of m was negative with very high correla¬ 
tion with G f . While AIC’s and BIC’s were high for GM C , they 
were lower for RM C and intermediate for MM C . Other diagnostic 
parameters were comparable. For L2, RM C did not converge to a 
solution in spite of various combinations of starting values. Pa¬ 
rameter estimates and diagnostic statistics of the other two mod¬ 
els indicated better characteristics in MM C . Trajectory of the 
fitted curves was more realistic in MM C (Fig. 3) while the resid¬ 
ual plots (Fig. 4) did not reveal any non-randomness in errors 
except in GM C . Bootstrap estimates of the parameters (Tables 6) 
were comparable with the respective model parameters. 



0 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 14 16 

Age (years) Age (years) 



10 15 20 25 30 35 40 45 50 55 10 15 20 25 30 35 40 45 50 55 


Predicted Predicted 

Fig. 2 Scatter plot of residuals against predicted values of the uncon¬ 
strained RM (A, B), GM (C, D) and MM (E, F) models fitted to data 
from two locations LI (A, C, E) and L2 (B, D, F) 



0 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 14 16 


Age (years) Age (years) 

Fig. 3 Scatter plot of girth of trees overlaid with function plots of 
constrained RM (A), GM (C, D) and MM (E, F) models fitted to data 
from two locations LI (A, C, E) and L2 (D, F). *RM model for this 
location did not converge to a solution and hence no overlaid function 
plot (Refer Table 5) 
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Table 4. Parameter estimates and goodness of fit statistics of the unconstrained models fitted to Hevea girth data from two locations 


Location 


Gf(cm) 

Estimate 

ASE 

95% ACI 

(L) 

95% ACI 

(U) 

G 0 (cm) 

Estimate 

ASE 

95% ACI 

(L) 

95% ACI 

(U) 

k Estimate 

ASE 

95% 

ACI (L) 

95% ACI 

(U) 

m Esti¬ 
mate 


RMu 

64.86 

2.05 

60.85 

68.88 

47.74 

0.09 

47.54 

47.93 

0.31 

0.03 

0.25 

0.38 

0.73 

Location 

GMu 

61.92 

0.6 

60.74 

63.1 

47.8 

0.09 

47.62 

47.98 

0.38 

0.01 

0.36 

0.39 

— 

(Ll) a 

MMu 

84.37 

2.14 

80.17 

88.57 

47.56 

0.09 

47.38 

47.73 

0.14 

0.01 

0.13 

0.15 

— 


RMu 

61.82 

3.06 

55.82 

67.81 

37.14 

0.09 

36.97 

37.32 

0.13 

0.03 

0.07 

0.18 

-0.85 

Location 

GMu 

52.56 

0.34 

51.9 

53.22 

37.12 

0.09 

36.94 

37.31 

0.34 

0.01 

0.32 

0.35 

— 

(L2) b 

MMu 

55.69 

0.52 

54.66 

56.71 

37.16 

0.09 

36.98 

37.34 

0.22 

0.01 

0.21 

0.24 

— 




95% 


Correlation 





Fit Statistics 




Location 



ACI 

95% ACI 

between Gf 












ASE 

(L) 

(U) 

& G 0 

G f & £ 

Gf&m 


RMSE 


LL 

AIC 

BIC 



RMu 

0.15 

0.45 

1.02 

-0.2 

-0.98 

-0.93 

RMSE 


LL 

AIC 

BIC 

MSC 


Location 

GMu 

— 

— 

— 

0.38 

-0.95 

— 

3.15 

0.94 

-6633.02 

13276.04 

13305.33 

2.84 


(Ll) a 

MMu 

— 

— 

— 

0.31 

-0.99 

— 

3.16 

0.94 

-6634.52 

13277.05 

13300.48 

2.84 



RMu 

0.26 

-1.35 

-0.35 

-0.32 

-0.99 

-0.94 

3.17 

0.94 

-6646.79 

13301.59 

13325.02 

2.83 


Location 

GMu 

— 

— 

— 

-0.64 

-0.93 

— 

3.92 

0.83 

-12558.02 

25126.03 

25158.1 

1.76 


(L2) b 

MMu 

— 

— 

— 

-0.69 

-0.97 

— 

3.94 

0.83 

-12584.51 

25177.03 

25202.68 

1.75 










3.92 

0.83 

-12564.33 

25136.65 

25162.31 

1.76 



ASE = asymptotic standard error; ACI = asymptotic confidence interval; L = lower; U = upper. RMSE = Residual mean square error; LL = Log likelihood; AIC = 
Akike information criterion; BIC = Bayesian information criterion; MSC = model selection criterion. a n = 34 trees, 76 cases tree" 1 ; b n= 48 trees, 94 cases tree" 1 


Table 5. Parameter estimates and goodness of fit statistics of the constrained models fitted to Hevea girth data from two locations. 


Parameter/Statistic 



Location (Ll) a 


Location (L2) b 



RM C 

GM C 

MM C 

RM C 

GM C 

MM C 


Gf (cm) 

Estimate 

148.69 

151.07 

138.84 


77.93 

93.18 


ASE 

1.73 

1.41 

0.53 


0.42 

0.40 


95% ACI (L) 

145.30 

148.30 

137.79 


77.11 

92.40 


95% ACI (U) 

152.09 

153.84 

139.88 


78.75 

93.96 

G 0 (cm) 

Estimate 

47.51 

47.53 

47.72 


35.50 

35.70 


ASE 

0.09 

0.11 

0.09 


0.07 

0.07 


95% ACI (L) 

47.33 

47.31 

47.54 


35.36 

35.57 


95% ACI (U) 

47.69 

47.74 

47.89 


35.64 

35.83 

k c 

Estimate 

0.04 

0.13 

0.07 


0.13 

0.07 

m 

Estimate 

-0.36 

__ 

__ 

The model did not 


__ 


ASE 

0.03 

_ 

_ 

converge to a solution 

_ 

_ 


95% ACI (L) 

-0.43 

__ 

__ 

in spite of various 

__ 

__ 


95% ACI (U) 

-0.29 

— 

— 

combinations of start- 

— 

— 

Correlation between Gf & G, 

0 

-0.09 

0.63 

0.79 

mg values 

-0.36 

-0.19 

Gf & m 


-0.96 

— 

— 


— 

— 

Fit Statistics 








RMSE 


3.20 

3.84 

3.26 


4.23 

4.13 

D 2 


0.94 

0.91 

0.94 


0.80 

0.81 

LL 


-6667.86 

-7144.24 

-6721.75 


-12906.43 

-12802.77 

AIC 


13343.73 

14294.48 

13449.49 


25818.86 

25611.55 

BIC 


13367.16 

14312.05 

13467.07 


25838.10 

25630.79 

MSC 


2.82 

2.45 

2.77 


1.61 

1.66 


ASE = asymptotic standard error; ACI = asymptotic confidence interval; L = lower; U = upper. RMSE = Residual mean square error; LL = Log likelihood; AIC = 
Akike information criterion; BIC = Bayesian information criterion; MSC = model selection criterion. a n = 34 trees,76 cases tree" 1 ; b n = 48 trees, 94 cases tree" 1 . 
c Baek calculated value based on the parameter’s relationship with duration of growth (T) as k = (2m+2)/30 for Richards, k = 4/30 for Gompertz and k = 2/30 for 
Monomolecular. 
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Table 6. Bootstrap estimates of model parameters of the unconstrained and constrained functions fitted to Hevea girth data from two locations 


Unconstrained model 


Location 


Gf (cm) 

Estimate 

95% BCI 

(L) 

95% BCI 

(U) 

G 0 (cm) 

Estimate 

95% BCI 

(L) 

95% BCI 

(U) 

k Estimate 

95% BCI 

(L) 

95% BCI 

m Estimate 
(U) 

95% BCI 95% BCI 

(L) (U) 


RM 

64.88 

61.5 

69.41 

47.73 

47.54 

47.93 

0.31 

0.25 

0.38 

0.73 

0.47 1.03 

Location 













fLlT 

GM 

61.94 

60.85 

63.14 

47.8 

47.62 

47.98 

0.38 

0.36 

0.39 

— 

— 


MM 

84.41 

80.59 

88.75 

47.56 

47.38 

47.74 

0.14 

0.13 

0.15 

— 

— 


RM 

61.7 

57.3 

70.85 

37.14 

37 

37.31 

0.13 

0.07 

0.19 

-0.84 

-1.31 -0.34 

Location 














GM 

52.57 

51.92 

53.23 

37.13 

36.95 

37.3 

0.34 

0.32 

0.35 

— 

-- 

(L2) a 














MM 

55.7 

54.74 

56.76 

37.16 

37 

37.35 

0.22 

0.21 

0.24 

— 

— 








Constrained model 





Location 


Gf (cm) 

95% BCI 

95% BCI 

G 0 (cm) 

95% BCI 

95% BCI 


95% BCI 

95% BCI 











m Estimate 







Estimate 

(L) 

(U) 

Estimate 

(L) 

(U) 


(L) 

(U) 




RM 

148.71 

145.55 

152.3 

47.51 

47.33 

47.69 

-0.36 

-0.42 

-0.3 



Location 













nn a 

GM 

151.07 

148.33 

153.78 

47.52 

47.32 

47.75 

— 

— 

— 




MM 

138.84 

137.82 

139.91 

47.72 

47.54 

47.9 

— 

— 

— 



Location 

RM 













GM 

77.92 

77.1 

78.78 

35.5 

35.36 

35.65 

-- 

-- 

— 



(L2) a 














MM 

93.17 

92.39 

93.94 

35.7 

35.57 

35.82 

— 

— 

— 




BCI = bootstrap confidence interval (with 2000 iterations); L = lower; U = upper. a n = 34 trees and 76 cases tree’ 1 in LI and 48 trees and 94 cases tree’ 1 in L2. 



10 20 30 40 50 60 10 20 30 40 50 60 

Predicted Predicted 


Fig. 4 Scatter plot of residuals against predicted values of the con¬ 
strained RM (A), GM (C, D) and MM (E, F) models fitted to data 
from two locations LI (A, C, E) and L2 (D, F) [*RM model for this 
location did not converge to a solution and hence no plot of residuals 
(Refer Table 5)] 

Mixed models 

In the mixed modelling approach (Tables 7 and 8), inclusion of 
random effects in both G f and G 0 resulted in significant reduction 
in RMSE, but did not result in removing the correlation between 


them. Addition of diagonal positive definite matrix to account for 
the variance-covariance structure of the random terms not only 
resulted in significant reduction in RMSE but also removed the 
correlation between G f and G 0 . The results were similar for both 
locations. Further, addition of variance and correlation structures 
did not result in significant improvement in model fittings. The 
residual plots of the models (Fig. 5) did not reveal any non ran¬ 
dom (systematic) distribution of errors but a seasonal component 
was evident. Growth profiles of the trees predicted for a growth 
period of 30 years (Fig. 6) using constrained random coefficient 
model (mixed model) indicated the superiority of LI for the stud¬ 
ied clone. Comparison of growth trajectory profiles of trees from 
both locations indicated that trees at L2 would take two to two 
and a half years more to attain the standard girth necessary for a 
decision to tap trees. 

6 
4 
2 
0 
-2 
4 


g -6 



2 

0 

-2 

-4 

-6 


Fig. 5 Scatter plot of residuals against predicted values of the con¬ 
strained random coefficient MM model fitted to data from two loca¬ 
tions LI (A) and L2 (B). 
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Table 7. Comparison of mixed-effects model performance with tree-level random effects components and different variance-covariance struc¬ 
tures for the constrained monomolecular model (MM C ) fitted to Hevea girth data for location 1. 


Specification 

Model 

df 

AIC 

LL 

RMSE 

r (G f vs. G 0 ) 

P LRT 

All parameters fixed 

Ml 

3 

13449.49 

-6721.75 

3.26 

0.79 

— 

Random effect in Gf 

M2 

4 

12916.3 

-6454.15 

2.88 

— 

— 

M2 vs. Ml 

— 

— 

— 

— 

— 

— 

0.001 

Random effects in Gf & G 0 

M3 

6 

9166.98 

-4577.49 

1.34 

0.914 

— 

M2 vs. M3 

— 

— 

— 

— 

— 

— 

0.001 

Variance -Covariance for random terms 

M3 + Diagonal positive-definite matrix 

M4 

5 

9221.27 

-4605.63 

1.34 

0.005 

— 

M3 vs. M4 

— 

— 

— 

— 

— 

— 

0.001 

M3 + General positive-definite matrix 

M5 

6 

9614.72 

-4576.36 

1.34 

0.914 

— 

M4 vs. M5 

— 

— 

— 

— 

— 

— 

0.001 

Variance structure for residuals 

M4 + Power of a variance covariate 

M6 

6 

9222.92 

-4605.46 

1.22 

0.006 

— 

M4 + Exponential of a variance covariate 

M7 

6 

9164.72 

-4576.36 

1.34 

0.912 

— 

M6 vs. M7 

— 

— 

— 

— 

— 

— 

ns 

Correlation structure for residuals 

M4 + Autoregressive process of order 1 (AR1) 

M4 + Continuous autoregressive process (AR1) for 

M8 

6 

4229.6 

-2108.8 

3.32 

0.649 

— 

a continuous time covariate (CAR1) 

M9 

6 

4229.6 

-2108.8 

3.32 

0.649 

— 

df = degrees of freedom; AIC = Akike information criterion; LL 

= Log likelihood; RMSE 

= Residual mean square error; r 

= correlation coefficient; 

P LRT = 

probability of likelihood ratio test, ns = not significant. 









Table 8. Comparison of mixed-effects model performance with tree-level random effects components and different variance-covariance struc¬ 
tures for the constrained monomolecular model (MM C ) fitted to Hevea girth data for location 2. 

Specification 

Model 

df 

AIC 

LL 

RMSE 

r (G f vs. Go) 

P LRT 

All parameters fixed 

Ml 

3 

25611.55 

-12802.77 

4.13 

-0.19 

— 

Random effect in Gf 

M2 

4 

24587.85 

-12289.93 

3.62 

-0.036 

— 

M2 vs. Ml 

— 

— 

— 

— 

— 

— 

0.001 

Random effects in Gf & G 0 

M3 

6 

17828.09 

-8908.55 

1.66 

0.786 

— 

M2 vs. M3 

— 

— 

— 

— 

— 

— 

0.001 

Variance -Covariance for random terms 

M3 + Diagonal positive-definite matrix 

M4 

5 

17873.27 

-8931.64 

1.66 

-0.001 

— 

M3 vs. M4 

— 

— 

— 

— 

— 

— 

0.001 

M3 + General positive-definite matrix 

M5 

6 

17829.09 

-8908.55 

1.66 

0.786 

— 

M4 vs. M5 

— 

— 

— 

— 

— 

— 

0.001 

Variance structure for residuals 

M4 + Power of a variance covariate 

M6 

6 

16898.34 

-8443.17 

16.77 

-0.006 

— 

M4 + Exponential of a variance covariate 

M7 

6 

16898.34 

-8443.17 

16.77 

-0.006 

— 

M6 vs. M7 

— 

— 

— 

— 

— 

— 

ns 

Correlation structure for residuals 

M4 + Autoregressive process of order 1 (AR1) 

M8 

6 

6379.97 

-3183.98 

4.09 

0.329 

— 

M4 + Continuous autoregressive process (AR1) for 








a continuous time covariate (CAR1) 

M9 

6 

6379.97 

-3183.98 

4.09 

0.329 

— 


df = degrees of freedom; AIC = Akike information criterion; LL = Log likelihood; RMSE = Residual mean square error; r = correlation coefficient; P |LRT 
probability of likelihood ratio test, ns = not significant. 
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Fig. 6 Growth profiles of Hevea trees predicted for a growth period 
of 30 years using the constrained random coefficient MM model for 
the two locations LI (A) and L2 (B). 

Discussion 

Fixed effects models 

Unconstrained functions underestimated the asymptotic maxi¬ 
mum that did not reflect the carrying capacity of the locations. 
This was expected because close to reality estimations are possi¬ 
ble only from a set of size-age measurements at every age until 
the attainment of plateau at maturity. The underestimations were 
due to partial sets of measurements made in the early growth 
phase of the trees. It has been observed at LI that untapped trees 
of clone RRII 105 attain 130 to 140 cm of girth after about 30 
years of growth (J. G. Marattukalam, (personal communication, 
2010). Underestimations can easily be attributed to the partial set 
of measurements made in early growth phase of the trees. Before 

settling on constraining (T), the option of constraining G f was 

considered and was discounted as inappropriate as it would not 
bring out productivity differences between the study sites. Con¬ 
straining duration of growth was the best option because it was 
possible to approximate a bench mark value for obtaining maxi¬ 
mum potential growth that would not exceed the biologically 
reasonable limits. Further, it was assumed that given a range of 
environments where rubber is being grown, fixing a bench mark 
value for the duration of growth would be prudent in standardiz¬ 
ing it for future studies. In addition, holding the growth period 
constant would enable more meaningful comparison of all 
growth quantities within, between, and among clones and loca¬ 


tions. 

Results from both unconstrained and constrained models indi¬ 
cated superiority of RM and MM compared to GM. Though RM 
yielded a fit as good as did MM, a slight edge of RM over MM 
was expected as the curve generated by an equation with more 
variables will nearly always have better diagnostic statistics sim¬ 
ply because it has additional parameters. In RM estimate of pa¬ 
rameter m was inconsistent between unconstrained and con¬ 
strained models both in sign and magnitude. While in RMy its 
value was positive for LI, negative for L2, but in RM C it was 
negative for LI and the model did not converge for L2. The cor¬ 
relation with G f and m was also high indicating the redundancy 
of the parameter. Further, negative values of m indicate absence 
of a lower asymptote and the inflection point, which are charac¬ 
teristic of the MM curve (Richards 1959; Rose & Charles- 
Edwards 1981; Zhao-gang & Feng-ri 2003; Lei & Zhang 2004). 
It has been reported that the Richards function will not fit sets of 
data which show insufficient curvature towards an upper asymp¬ 
tote (Hunt 1982). Further, Fekedulegn et al. (1999) reported that 
investigations of the differential forms and second derivatives of 
the equation have indicated the unsuitability of the function to 
model data that does not encompass the entire range of the life 
cycle (i.e., juvenile, adolescent, mature and senescent stages of a 
biological variable). Some authors consider that the additional 
parameter used in RM has no obvious biological interpretation 
and it is numerically unstable and often poses problems in fitting 
the model (Thomley and Johnson 1990; Zeide 1993). 

Random effects models 

In the random coefficient models, both G f and G 0 appeared to be 
influenced by tree level effects. Inclusion of a diagonal positive 
definite matrix to account for the variance-covariance structure 
of the random terms removed the correlation between the ran¬ 
dom effects. The results were similar for both locations. No fur¬ 
ther improvement could be achieved by including variance mod¬ 
elling to account for heteroscedasticity and within subject corre¬ 
lation in the presence of variance-covariance structures for the 
random effects. Gregoire and Shabenberger (1996) pointed out 
that additional within subject correlation is often negligible in the 
presence of random subject effects. 

Overall assessment 

In plant growth modelling, the use of fitted curves is of very high 
importance because the fitted curves can provide the experi¬ 
menter with a clearer perception of the reality of plant growth 
when a series of observational data are disturbed by random 
errors. In most cases the curve is needed not only to provide a 
convenient re-description of the primary data but as an instru¬ 
ment for derivation of the quantities involved in plant growth 
analysis (Causton and Venus 1981; Hunt 1982; France and 
Thornley 1984). In Hevea , growth measured in terms of increase 
in girth of the main stem remains the most important variable of 
interest to commercial planters and managers. Girth at age seven 

^ Springer 





374 


Journal of Forestry Research (2012) 23(3): 3654375 


is the foremost factor taken into account for evaluating the pro¬ 
gress of growth and attainment of maturity for tapping. Girth and 
rate of girth increase are also used in experimental work to assess 
the growth performance of new planting materials and effects of 
cultural treatments on growth (Shorrocks et al. 1965; Chan- 
drashekar et al. 1998; Chandrasekhar et al. 2005; Chandrasekhar 
2007). 

From my results and discussion it is clear that MM is the pre¬ 
ferred candidate model for studying the girth-age relationships in 
Hevea trees. The asymptotic maximum obtained for this model 
was closer to the observed information for the clone at LI. Addi¬ 
tionally, values of allometric parameter m<0 and concave curva¬ 
ture of the scatter plots of observed data also indicate the suit¬ 
ability of MM. Goodness of fit, used alone, could lead to choos¬ 
ing an inappropriate growth function: The extent to which a 
growth function produces reasonable biological estimates must 
remain a primary factor in model selection (Vanclay and Skovs- 
gaard 1997; Caillet et al. 2006). This model, though not com¬ 
monly employed in forest growth studies, is widely used in mod¬ 
elling fast-growing forest trees (Amaro et al. 1998; Zhao-gang & 
Feng-ri 2003; Lei & Zhang 2004). Thus it appears that in Hevea 
trees, photo synthetic activity decreases as size [G(t)] increases 
and constructive metabolism exceeds destructive metabolism. As 
a result, the growth rate is maintained at maximum value at t 0 , 
then decreases until the final stage at dG/dt>0 resulting in a yield 
curve with no period of accelerating growth. The decelerating 
growth takes place in the beginning and lasts until the final limit¬ 
ing size, G f (Richards 1959; Zhao-gang & Feng-ri 2003; Lei & 
Zhang 2004). 

Conclusions 

The MM model is the most suitable function for study of the 
increase in girth of Hevea trees. Further, for subject specific 
modelling, inclusion of a diagonal positive definite matrix to 
account for the variance-covariance structure of the random 
terms would take care of both heteroscedasticity and within sub¬ 
ject correlation in the repeated measures. The model presented is 
a generic one ignoring the influences of various other factors. 
Testing of this model utilising data from contrasting climatic 
conditions and incorporation of modifier functions of weather 
factors into the growth equation should yield more information 
on the utility of the function to model Hevea growth. Further 
studies are in progress to incorporate seasonality of growth (both 
wet and dry seasons), clonal variation, site index (that reflect the 
nutrient availability and average climatic conditions), environ¬ 
mental influences (rainfall, temperature, relative humidity, soil 
factors, and others) and management practices (plating density, 
fertilisation regimes, and others) in the growth model. 

References 

Amaro A, Reed D, Tome M, Themido I. 1998. Modelling dominant height 

growth: Eucalyptus plantations in Portugal. Forest Science, 44: 3746. 

4^ Springer 


Baulkwill WJ. 1989. The history of natural rubber production. In: C.C. Web¬ 
ster and W.J. Baulkwill (eds), Rubber. UK: Longman Scientific and Tech¬ 
nical, pp. 1-56. 

Caillet GM, Smith WD, Mollet HF, Goldman KJ. 2006. Age and growth 
studies of chondrichthyan fishes: the need for consitency in terminology, 
verification, validation and growth function fitting. Environmental Biology 
of Fishes, 77: 211428. 

Causton DR, Venus JC. 1981. The biometry of plant growth. Edward Arnold, 
UK. 

Chandrasekhar TR, Alice J, Varghese YA, Saraswathyamma CK, Vijayaku- 
mar KR. 2005. Girth growth of rubber (Hevea brasiliensis ) trees during the 
immature phase. Journal of Tropical Forest Science , 17: 399415. 

Chandrashekar TR. 2007. Weather and growth performance of young rubber 
trees (Hevea brasiliensis). Ph. D. Thesis, Mahatma Gandhi University, 
Kerala, India. 

Chandrashekar TR, Nazeer MA, Marattukalam JG, Prakash GP, Anna- 
malainathan K, Thomas J. 1998. An analysis of growth and drought toler¬ 
ance in rubber during the immature phase in a dry subhumid climate. Ex¬ 
perimental Agriculture, 34: 287400. 

Delignette-Muller ML, Baty F. 2010. Use of package nlstools to help the fit 
assess the quality of fit of a gaussian nonlinear model. Available at 
http://cran.r-project.org/web/packages/nlstools. [accessed 27 December 
2010 ]. 

El-Shaarawi AH, Piegorsch WW. 2002. Encyclopedia of environmetrics 
(Vol.l), UK: Wiley, pp 32-41.. 

Fekedulegn D, Mac Siurtain MP, Colbert JJ. 1999. Parameter estimation of 
nonlinear growth models in forestry. Silva Fennica, 33: 327436. 

France J, Thornley JHM. 1984. Mathematical models in agriculture. UK: 
Butterworths, pp 223435. 

Gregoire TG, Schabenberger O. 1996. A non-linear mixed effects model to 
predict cumulative bole volume of standing trees. Journal of Applied Sta¬ 
tistics, 23: 257471. 

Gregoire TG, Brillinger DR, Diggle PJ, Russek-Cohen E, Warren WG, Wolf- 
inger RD. 1997. Modelling longitudinal and spatially correlated data. 
Springer-Verlag, New York. 

Hunt R. 1982. Plant growth curves, the functional approach to plant growth 
analysis. London: Edward Arnold. 

Lei YC, Zhang SY. 2004. Features and partial derivatives of Bertalanffy- 
Richards growth model in forestry. Non linear analysis: Modelling and 
control, 9: 6543. 

Liu Zhao-gang, Li Feng-ri. 2003. The generalized Chapman-Ricahrds func¬ 
tion and applications to tree and stand growth. Journal of Forestry Re¬ 
search, 14: 1646. 

Ogle DH. 2010. Von Bertalanffy growth model vignette. Available from 
http://www.ncfaculty.net/dogle/fishR/gnrlex/VonBertalaffy/VonBertalanffy 
.pdf. [Accessed 27 December 2010]. 

Philip MS. 1994. Measuring trees and forests. 2 nd Edition. Wallingford, UK: 
CAB International,. 

Pinheiro JC, Bates DM. 2000. Mixed effects models in S and S-PLUS. New 
York: Springer-Verlag. 

Potty SN, Kodandaraman R, Mathew M. 1980. Field upkeep. In: P.N. Rad- 
hakrishna Pillay (ed), Hand book of Natural Rubber Production in India. 
Rubber Research Institute of India, Kottayam, India, pp 135456. 

Punnoose KI, Kodandaraman R, Philip V, Jessy MD. 2000. Field upkeep and 
intercropping. In: P.J. George and C.K. Jacob (eds), Narural Rubber Agro¬ 
management and Crop Processing. Kottayam, India: The Rubber Research 




Journal of Forestry Research (2012) 23(3): 365^75 


375 


Institute of India, pp 149469. 

Pushpadas MV, Ahammed M. 1980. Nutritional requirements and manorial 
recommendations. In: P.N. Radhakrishna Pillay (ed), Hand book of Natural 
Rubber Production in India. Kottayam, India: Rubber Research Institute of 
India, pp 159485. 

R Development Core Team. 2010. R: A language and environment for statis¬ 
tical computing. R Foundation for Statistical Computing, Vienna, Austria. 
ISBN 3-900051-07-0. Available from http://www.r-proiect.org [accessed 
31 January 2011]. 

Ratkowsky DA. 1988. Handbook of Nonlinear Regression Models. Marcel 
Dekker, NewYork. 

Richards FJ. 1959. A flexible growth function for empirical use. Journal of 
Experimental Botany, 10: 290-300. 

Rose DA, Charles-Edwards DA. 1981. Mathematics and plant physiology. 
Academic press, London. . 

Salas C, Garcia O. 2006. Modelling height development of mature Not- 
hofagus obliqua. Forest Ecology and Management, 229: 1-6. 

Seber GFA, Wild CJ. 2003. Nonlinear Regression. USA: John Wiley & Sons 
Inc.,. 

Shifley SR, Brand GJ. 1984. Chapman-Richards growth function constrained 
for maximum tree size. Forest Science , 30: 10664070. 

Shorrocks VM, Templeton JK, Iyer GC. 1965. Mineral Nutrition, growth and 


nutrient cycle of Hevea brasiliensis III. The relationship between girth and 
shoot dry weight. Journal of Rubber Research Institute of Malaya, 19: 
83-92. 

Thornley JHM, France J. 2007. Mathematical Models in Agriculture, 2nd ed. 
Wallingford, UK: CAB International. 

Thornley JHM, Johnson IR. 1990. Plant and crop modelling: a mathematical 
approach to plant and crop physiology. Oxford, UK: Oxford University 
Press, 

Vanclay JK. 1994. Modelling forest growth and yield. Wallingford, UK: CAB 
International,. 

Vanclay JK, Skovsgaard JV. 1997. Evaluating forest growth models. Ecologi¬ 
cal Modelling, 98: 142. 

Vonesh EF, Chinchilli VM. 1997. Linear and nonlinear models for the analy¬ 
sis of repeated measurements. New York: Marcell Dekker. 

Webster CC, Paardekooper EC. 1989. The botany of the rubber tree. In: C.C. 
Webster and W.J. Baulkwill (eds), Rubber. England: Longman Scientific 
and Technical, pp 57-84. 

Winsor CP. 1932. The Gompertz curve as a growth curve. Proceedings of the 
Nattionla Academy of Sciences, 16: 1-8. 

Zeide B. 1993. Analysis of growth equation. Forest Science, 39: 594-616. 


4^ Springer 




